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Abstract: We present a subtraction method utilizing the A-jettiness observable, Tn, to 
perform QCD calculations for arbitrary processes at next-to-next-to-leading order (NNLO). 
Our method employs soft-collinear effective theory (SCET) to determine the IR singular 
contributions of A-jet cross sections for Tn —> 0, and uses these to construct suitable 
7iv-subtractions. The construction is systematic and economic, due to being based on 
a physical observable. The resulting NNLO calculation is fully differential and in a form 
directly suitable for combining with resummation and parton showers. We explain in detail 
the application to processes with an arbitrary number of massless partons at lepton and 
hadron colliders together with the required external inputs in the form of QCD amplitudes 
and lower-order calculations. We provide explicit expressions for the 7Ar-subtractions at 
NLO and NNLO. The required ingredients are fully known at NLO, and at NNLO for 
processes with two external QCD partons. The remaining NNLO ingredient for three or 
more external partons can be obtained numerically with existing NNLO techniques. As an 
example, we employ our results to obtain the NNLO rapidity spectrum for Drell-Yan and 
gluon-fusion Higgs production. We discuss aspects of numerical accuracy and convergence 
and the practical implementation. We also discuss and comment on possible extensions, 
such as more-differential subtractions, necessary steps for going to N 3 LO, and the treatment 
of massive quarks. 
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1 Introduction 


The precise knowledge of QCD corrections is a key ingredient for interpreting the data 
from collider experiments. In hadronic collisions, the inclusive QCD cross section for the 
production of a final state X can, if the hard scale Q associated with X is large enough, 
be obtained in terms of a perturbatively calculable partonic cross section convolved with 
parton distribution functions (PDFs). 

Perturbative calculations performed using the leading order (LO) term in a s typically 
suffer from large theoretical uncertainties due to missing higher-order perturbative correc¬ 
tions. Often, next-to-leading order (NLO) is the first order at which the normalization and 
in some cases the shape of cross sections can be considered reliable. As such, this level of 
accuracy has become standard for comparing with data from the LHC. For some processes 
the experimental uncertainties are becoming so small, or the perturbative uncertainties at 
NLO are still so large, that next-to-next-to-leading order (NNLO) computations are called 
for. 

For many important benchmark processes, the required virtual amplitudes are known 
at NNLO. However, as is well known, the computation of the full cross sections beyond 
leading order is complicated by infrared (IR) divergences - explicit divergences in virtual 
amplitudes, and divergences in the phase-space integration over the real-emission ampli¬ 
tudes in regions where particles become soft or collinear to other particles. These diver¬ 
gences only cancel after integrating the real-emission amplitudes over the phase space of 
unresolved particles and adding the result to the virtual loop amplitudes order by order. 

To handle these divergences in practice one typically makes use of some subtraction 
method. That is, one subtracts terms from the real emission contributions that reproduce 
the IR soft and collinear behaviour of the real emissions, which then allows the phase-space 
integral of the full amplitude minus the subtraction terms to be performed numerically in 
d = 4 dimensions, giving a finite result. The subtracted terms have to be sufficiently 
simple that they can be integrated over the phase space of emitted particles in d = 4 — 2e 
dimensions. They are then added back to the virtual contributions, where they cancel the 
explicit l/e n IR poles. 

The goal of typical NLO subtraction schemes like FKS subtractions [1-3] or CS sub¬ 
tractions [4-6] is to construct subtraction terms that reproduce the correct IR-singular 
behaviour of the full real-emission amplitude point-by-point in phase space. Over the past 
decade enormous effort has been devoted to extend such local subtraction methods to 
NNLO using different approaches [7-38]. This extension is very involved due to the many 
overlapping singularities at NNLO, which have to be isolated by appropriate phase-space 
parameterizations. At the same time, the subtractions have to remain simple enough that 
the l/e n IR poles can be extracted from the integrated subtractions. 

The basic idea of our method, which we call AQjettiness subtractions, is to use a 
physical jet-resolution variable Tn to control the infrared behaviour of the cross section. 
The key point is that, if the (factorized) structure of the leading contribution to the Tn- 
differential cross section in the IR limit Tn — > 0 is known, the singular part can often be 
determined analytically and used to construct an IR subtraction term. A major advantage 
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of using a physical observable is that the differential and integrated subtraction terms 
are then equivalent to the singular limits of a physical cross section, which can indeed be 
significantly easier to calculate than the full cross section. A well-known example of such 
a physical subtraction scheme is the (^-subtraction method for color-singlet production in 
hadron collisions [39], which has been successfully applied to a variety of processes [40-47]. 
(It has also been suggested that this method can be applied to compute heavy-quark pair 
production at NNLO [48, 49].) Our A-jettiness subtraction method generalizes this to 
arbitrary numbers of QCD partons in the initial and final state. It employs the IV-jettiness 
global event shape [50] as the physical A-jet resolution variable. In this paper, we limit 
ourselves to massless quarks; the extension to massive quarks is in principle possible and 
commented on in section 5. 

The key feature of A-jettiness is that it has very simple factorization properties in the 
singular limit. The factorization theorem for the A-jettiness cross section is known [50— 
52] from soft-collinear effective theory (SCET) [53-58]. It can be used to systematically 
compute the leading singular contributions (thus determining the subtraction terms) by 
performing standard fixed-order calculations of soft and collinear matrix elements in SCET. 
At NLO, all necessary ingredients have been known for some time, and by now, essentially 
all necessary NNLO ingredients are available. For processes with hadronic initial states a 
key ingredient that has become available recently are the two-loop quark and gluon beam 
functions [59, 60]. 

The price one has to pay for using a single physical observable to describe the IR is that 
the subtraction does not act point-by-point in phase space, but only on a more global level 
after a certain amount of phase-space integration has been carried out. In essence, the large 
number of terms in a fully local subtraction method are projected onto a single, nonlocal 
subtraction term. In practice, this means that the numerical convergence may be slower 
than for the fully local case. However, this is compensated by the significant reduction 
in complexity of the subtractions. Furthermore, as we will discuss, it is possible to make 
the subtractions step-by-step more local by making the A-jettiness cross section more 
differential in additional variables. This is again possible by using SCET to factorize and 
calculate the singular contributions of more differential cross sections (see e.g. refs. [52, 61— 
65]). 

There are several important benefits of using a physical observable as jet resolution 
variable, as already emphasized in refs. [66]. It allows one to directly reuse the existing NLO 
calculations for the corresponding A + 1-jet cross sections, and the resulting NNLO calcu¬ 
lation is automatically fully-differential in the Born phase space. Moreover, the calculation 
will be in a form which makes it directly suitable to be combined with higher-order resum¬ 
mation as well as parton showers by using the general methods developed in refs. [66, 67]. 

The idea of using A-jettiness as an A-jet resolution variable is not new. In fact, this 
is what largely motivated its invention in the first place. It is already utilized in essentially 
the same context as here in the Geneva Monte-Carlo program [67] . For color-singlet pro¬ 
duction, the A-jettiness subtraction method reduces to an analogue of qx subtractions [39] 
with an alternative physical resolution variable. The differential version as a subtraction 
was used at NLO in ref. [68]. 
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In its simplest form as a phase-space slicing, the IV-jettiness subtraction method has 
been successfully applied already to calculate the top quark decay rate at NNLO [69]. 1 
While this work was being finalized, this method was also suggested and applied to the 
NNLO calculations of pp W/H + jet in Refs. [72, 73]. These results clearly highlight 
the usefulness of the slicing method, even for complex 2 —> 2 processes with three colored 
partons. 

In this work we give a general description of how IV-jet resolution variables, and specif¬ 
ically IV-jettiness, can be used as subtraction terms to compute fixed-order cross sections. 
In section 2, we discuss how the IR singularities in QCD cross sections are encapsulated 
by an IV-jet resolution variable. We demonstrate that this naturally leads to subtraction 
terms for fixed-order calculations, and show how these can be used in phase-space slicing, 
as done in refs. [69-73], and as differential subtractions, generalizing (^-subtractions [39]. 
In section 3, we review the definition of IV-jettiness and its general factorization theorem 
for IV-jet production. We show how the subtraction terms are defined in terms of func¬ 
tions in the factorization theorem. We explicitly construct the subtraction terms at NLO 
and NNLO for generic IV-parton processes. We also discuss the extension to N 3 LO and 
to more-differential subtractions. In section 4, we discuss how these subtractions may be 
implemented in parton-level Monte-Carlo programs. We also show results for Drell-Yan 
and gluon-fusion Higgs production at NNLO and use these as an example to discuss some 
of the numerical aspects. We conclude in section 5. 

2 General Formalism 
2.1 Notation 

We denote the IV-jet cross section that we want to compute by cr(X). Here, X collectively 
stands for all differential measurements and kinematic cuts applied at Born level. In 
particular, it contains the definitions of the N identified signal jets in cr(X) and all cuts 
required to stay away from any IR-singularities in the IV-parton Born phase space. 

The cross section at leading order (LO) in perturbation theory can then be written as 

a LO (X) = Jd<S> N B N ($ N )X(<f> N ), (2.1) 

where the measurement function X($at) implements X on an IV-parton final state. The 
Born contribution, is given by the square of the lowest-order amplitude, M®, for 

the process we are interested in, 2 

B n (<S> n ) = J2 \^\$n)\ 2 or B N (Q N ) = f a f b ( 2 - 2 ) 

color color 

where Tat denotes the complete dependence of the amplitude on the external state (includ¬ 
ing all dependence on momentum, spin, and partonic channel). For hadronic collisions, 

X A similar slicing method utilizing heavy-quark effective theory was also used in ref. [70, 71] to perform 
the fully-differential NNLO calculation for e + e“ —> tt. 

2 For a tree-level process, is given by the sum of the relevant tree-level diagrams. For a loop-induced 

process, like gg — > H, it is the sum of the relevant lowest-order IR-finite loop diagrams. 
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the PDFs f a ,b are included in Bn{&n) and 4 >tv also includes the corresponding momentum 
fractions x a ,b- Correspondingly, the integral over d^Ar in eq. (2.1) includes all phase-space 
integrals and sums over helicities and partonic channels. For simplicity, we also absorb into 
it flux, symmetry, and color and spin averaging factors. We use N to denote the number 
of strongly-interacting partons in the final state. There can also be a number of additional 
nonstrongly interacting final states at Born level, which are included in 4 ?tv but we suppress 
for simplicity. 


2.2 Singular and nonsingular contributions 


Any IV-jet cross section cr(X) can also be measured differential in a generic IV-jet resolution 
variable Tn, which we write as da(X)/dl~N ■ Then cr(X) may be written as 


a(X) 



da(X) 
d Tn 



da{X) 
d T n 


+ 



dcr(X) 
d Tn 


(2.3) 


dividing the more differential cross section into the region 0 < Tn < T§ ut and the re¬ 
gion Tn > T^ ut . For Tn to be an IV-jet resolution variable it must satisfy the following 
conditions: 


Tn($n) = 0, Tn(®>n+i)>0, Tn{$>n+i —> d>iv) —> 0. (2.4) 


In words, Tn must be a physical IR-safe observable that resolves all additional IR-divergent 
real emissions, such that the cross section d<r(X)/d7iv is physical and IR finite for any 
Tn > 0, and the IR singular limit corresponds to Tn -A 0. 3 Hence, we have 


d<r LO (X) 


d Tn 


= a 


LO 


(■ X)6(T n ), 


= 0(a s 


(2.5) 


Tn> 0 


1 dcr(A') 

<r LO (X) dT N 

We use the convention that Tn is normalized to be a dimension-one quantity, and for 
convenience we also define the dimensionless quantities 


= Tn 
' Q 


*-rent 

cut _ ' N 

Q 


( 2 . 6 ) 


Here, Q is a typical hard-interaction scale of the Born process (whose precise choice however 
is unimportant). For example, canonical choices would be Q = E cm for e + e _ —> jets, 
Q = for Drell-Yan pp —>■ V —>■££, Q = mu for gg —> H , and Q = for pp —> dijets. 

We define the “singular” part of the Tn spectrum to contain all contributions that are 
singular in the Tn —> 0 limit, i.e., all contributions which are either proportional to 6(Tn) 
or that behave as ln n (r)/r for r —> 0. It can be written as 


daS 1 J (X) = C-i(X) 5(t) + Cn(X) £n(r) , (2.7) 

T n> 0 

3 For particular definitions of Tn, there could also be regions of 4>>iv+i (far) away from any IR singu¬ 
larities where Tn is small or vanishing. Such regions do not pose a problem and are irrelevant for our 
discussion. The typical example for Tn = qr at NNLO are contributions from two hard real emissions that 
are back-to-back such that qr —> 0. Another generic example are regions where two partons are collinear 
that cannot arise from a QCD singular splitting. Such cases can be avoided by defining Tn in a flavor-aware 
way. 
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where the C n (r) are the usual plus distributions. For a suitable test function /(r): 


C-nij) = 


dr £, 


8(t) ln n (r) 


ln n (r) 


ln n+1 (r 


i(T)f(r)= [ dT — v ' ' [/(t) - /(0)] + /(0) 

Jo T n- 1 - i 


cut^ 


( 2 . 8 ) 


This logarithmic structure of the singular contributions directly follows from the IR singular 
structure of QCD amplitudes, the KLN theorem, and the fact that 7 at is an IR-safe physical 
observable. Since the infrared limit of the QCD amplitudes, and hence the IR singularities, 
depends only on the lower-order phase space, the singular coefficients C n only depend on 
the underlying 4 >at. That is, 

= (2.9, 

We can therefore consider the singular distributions directly as a function of the full 
and independently of the specific measurement X, 


der sing (^>jv) 

dr 


C-i(<5jv) <5(r) + ^ C n (& N ) C n (r) 

n> o 


m >0 L 


2m— 1 




n =0 



( 2 . 10 ) 


In the second line, we have expanded the singular coefficients in a s . At LO, the only 
nonzero coefficient is 

C^ n ) = B n (^ n ), (2.11) 

so at LO the singular spectrum reproduces the LO cross section, consistent with eq. (2.5), 

= ci°j (X) 5(T n ) = a LO (X) 5{T n ) . (2.12) 

d In 

At NLO, the coefficients C_i j o,i( < l ) Ar) are nonzero, while at NNLO, the coefficients 
C—i,o, 1 , 2 , 3 (^iv) contribute. 

Writing the singular spectrum in terms of plus distributions as in eqs. (2.7) and (2.10) 
precisely encodes the cancellation between real and virtual IR divergences. The C_i coef¬ 
ficient contains the finite remnant of the virtual contributions after the real-virtual cancel¬ 
lation has taken place. By itself, it is not unique, but depends on the boundary conditions 
adopted in the definition of the plus distributions, which is encoded in the choice of r (the 
choice of Q ). Changing the boundary conditions is equivalent to rescaling the arguments 
of the plus distributions according to (see e.g. ref. [74]) 

” /rA ln n+1 A 

XC n (Xr) = J2 [ k )^ k X£n-k(r) + —- 5(t) . (2.13) 

k =0 ^ ' U + 
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While this rescaling moves contributions between different C n , it does not change the 
overall 1/Tn scaling, which implies that the sum of all terms in eq. (2.7) is unique 4 5 and 
in fact independent of the choice of Q.~ J Once the singular spectrum is written in terms 
of distributions as in eq. (2.7), one can easily integrate it up to Tn < T^ ut to obtain the 
singular cumulative distribution (or cumulant in short) 


0-sing 



dosing ( X ) 

d T n 


C- 1 (X) + J2 C n(X) 

n> 0 


ln n+1 (r cut ) 

n + 1 


(2.14) 


The “nonsingular” contributions are defined as the difference between total and sin¬ 
gular contributions, 

da nons {X) _ d ,<r(X) da sin s(X) 

d7 n dTv d7)v ’ 

r t n' 1 * j^.nons/ v\ 

a nons (X, T§ ut ) = / dT N ’ = a{X, T^) - a sing (X, T^) . (2.15) 

jo d / at 

They start at 0(a s ) relative to cr LO (A') (which is part of dcr smg ). By definition of the 
singular terms, the nonsingular spectrum contains at most integrable singularities for Tn — > 
0, the largest terms being da nons (X)/dTN ~ a”ln 2n (r). Equivalently, the nonsingular 
cumulant behaves for T^ ut —> 0 as 

cr nons (A:, 7/v ut ->■ 0) ~ r cut a? In 2n ( r cut ) -> 0 . (2.16) 


Hence, also the underlying matrix-element contributions yielding the nonsingular terms 
can be safely integrated in the infrared. 


2.3 7yv-subtractions 

Up to this point, the decomposition of a cross section into singular and nonsingular terms is 
just notation and holds for any Tn- The key point of the 7)v-subtraction method is that if 
we have analytic control of the singular Tn dependence, we can turn the singular spectrum 
dcr smg (Ai)/d7iv and its integral cr sing (A", T^ ut ) into subtractions, as discussed next. This 
requires that for some IV-jet resolution variable Tn, the underlying coefficients C u (<&n) in 
eq. (2.10) can be determined explicitly/’ In particular, the ability to explicitly compute 
C_i (Ttv) is precisely equivalent to being able to compute the integrated subtractions in a 
classical subtraction method. All these conditions are satisfied for IV-jettiness, as we will 
discuss in section 3. 

4 It is unique in the sense that it has the minimal Tn dependence, only containing In 71 (Tn ) /Tn ■ One 
could in principle include some subleading Tn dependence in the coefficients, if this turns out to be useful 
or convenient. This would move some contributions between the singular contributions and the nonsingular 
remainder in eq. (2.15). 

5 The actual physical scales appearing together with Tn in the logarithms are set by the hard Born 
kinematics. The reason to think of Q as a typical hard scale is that this provides the natural power 
suppression of the nonsingular terms. 

b They do not necessarily have to be known fully analytically, and in general they will not be. All we 
really need is a sufficiently fast way to compute their numerical values for given 4 >jv to in principle any 
desired accuracy. 
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2.3.1 7/v-slicing 

If the singular contributions for a given Tn are known, we can use T^ ut to divide the phase 
space into two regions: Tv < T§ ut and Tn > 7)!) ut . Taking T^ ut —> Ts = 5ir Q, where 
<5 ir = Ts/Q is an (in-principle) arbitrarily small IR cutoff, the singular terms will numer¬ 
ically dominate the nonsingular for Tn < Tjy ut - In fact, since the nonsingular cumulant 
a nons (X,Ts) is of 0(Ts/Q ) = 0(<5ir), we can neglect it in this limit. Hence, we get 





dcj(X) 

d Tn 


+ 



da(X) 
d Tn 


= <r sing (X, Ts) + [ d T n 
■JT S 


da(X) 
d Tn 


+ 0(8 m ). 


(2.17) 


This is precisely a phase-space slicing method, which we will call Tv-slicing. Calculating 
&(X) to N ra LO in this way requires determining a smg (X,Ts) to N n LO, which includes the 
N n LO virtual contributions. Beyond that, since the Tn spectrum only starts at 0(a s ) 
relative to cr(X), the problem is reduced to the N n-1 LO calculation for the cross section 
da(X)/dT\r for Tn > Ts- Furthermore, if an N n_1 LO calculation is available, the slicing 
only needs to be performed for the pure N n LO terms. 


2.3.2 Differential Tv-subtractions 

It is instructive to rewrite the Tv-slicing in eq. (2.17) in the form of a subtraction as follows, 


a(X) = a si ^(X,T oS ) + 



d a(X) 



da sins (X) 

dT N 


+ 0(S m ). 


(2.18) 


This reorganization shows that the integral of the singular spectrum acts as a global sub¬ 
traction for the integrated full spectrum, while the cumulant <7 sing (X, T 0 s) is the corre¬ 
sponding contribution of the virtual terms (sitting at Tn = 0) plus the integrated subtrac¬ 
tion. The value of T 0 g is arbitrary and exactly cancels between the first and third terms. It 
determines the upper limit in Tn up to which the subtractions are used. The subtraction 
term in this case is maximally nonlocal, as it is applied after all phase-space integrations. 
Hence, one would naively expect the numerical cancellations to be maximally bad. This 
also shows that Ts really is an IR cutoff below which only the singular (subtraction) terms 
are used, due to limited numerical precision. 

Looking at eq. (2.18), we can also move the singular spectrum underneath the Tn 
integration, 


dT/v 


<J(X) = a sm ^X,T oS ) + / d T n 

■JTs 

rTft j^-nons 

= a si ^(X,T oS )+ / dT N — 

JTs d 'N 


da(X) d a sing (X) 


dT N 


0(T n < Toff) 


+ 0(5i R ) 


PO 


+ 


'Ts 


da(X) 

dTv 


+ OPir) • 


(2.19) 


which turns the singular spectrum into an actual subtraction which is local (point-by¬ 
point) in Tn- It is of course still nonlocal in the remaining real radiation phase space. To 
use eq. (2.19), one now has to explicitly calculate the singular differential spectrum. This 
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requires essentially no additional effort, since the required singular coefficients are the same 
as in (T sin §(X,7^ ut ). 

Writing it as in the second line of eq. (2.19) shows explicitly that the numerical integral 
over Tn now only encounters an integrable singularity for 7jv —> 0 since the integrand is 
precisely the nonsingular contribution. This turns Ts into a purely technical cutoff for the 
numerical integration, which is only necessary because the integrand is still given by the 
difference of two diverging integrands. Finally, we note that the neglected contributions 
due to the numerical IR cutoff Ts are precisely the same as in eq. (2.17) for the same value 
of Ts- The numerical error introduced by such a cutoff is discussed in the next section. 

We stress that a technical IR cutoff analogous to <5 ir exists in any numerical fixed- 
order calculation using subtractions, since the QCD amplitudes (and their subtractions) 
become arbitrarily large in the IR. Below the cutoff, the full QCD amplitudes are always 
approximated by the subtraction terms, so that below the cutoff only the integral of the 
subtraction is used, while the nonsingular cross section below the cutoff is power suppressed 
by e>iR and neglected. 

Finally, note that separating the spectrum or cumulant into its singular and nonsin¬ 
gular parts, as we have done here, is in fact very well known and routinely used when 
performing the higher-order resummation for an IR-sensitive observable Tn- In this con¬ 
text, the singular contributions are resummed to all orders in a s and a given logarithmic 
order, while eq. (2.15) is used to determine the nonsingular contributions. At NNLO, this 
utilizes the result for da(X)/dT\r obtained from the NLO N + 1-jet calculation and the 
NNLO singular contributions obtained from the NNLL' resummation of Tn- In section 3 we 
will employ the same techniques to compute directly the fixed-order singular contributions 
without resummation. This also makes it clear that if desired any NNLO calculation per¬ 
formed in this way can be straightforwardly improved with the corresponding higher-order 
resummation in Tn- 


2.3.3 Estimating numerical accuracy 

We can judge the numerical accuracy of the 7Ar-slicing and differential 77v-subtractions 
using some simple scaling arguments. First, it is important to quantify the effect of the IR 
cutoff diR. Using IV-jettiness as an example, at N n LO relative to the Born cross section, the 
most dominant singular terms in the spectrum and the cumulant are, for a given partonic 
channel, 


^- LO ES(g)”(-Earo)N 2 „_ lW + .. 

n> 1 i 


o(V?')=cr LO '£ 


1 /a.,\ n 


n> 1 


n\ 


! \47t 


Y, C i Fo) ln 2 n (r cut ) + 


( 2 . 20 ) 


Here, To = 4 is the one-loop coefficient of the cusp anomalous dimension, C t = Cp for 
quarks and Ci = Ca for gluons, and the ellipsis denote terms with fewer powers of loga- 
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rithms at each order in a s .' Correspondingly, the leading nonsingular term in the cumulant 
has the form 


<7 


n> 1 


(7S ut ) = * L0 E h (S)'"ca, (■- E 1 c<r„ 


47T/ 


cut 


ln 2n_1 (r cut ) 


+ ■■■ 


( 2 . 21 ) 


The coefficient Cnols is not known in general, but we take cions = 1 here, which is the 
correct value for 2-jettiness in e + e _ (i.e. thrust). 

We denote the missing nonsingular contribution due to approximating the full result 
by the singular contributions below 7 jv < 7$ by Actir(<5ir) and expand it in a s as 

u nons (Ts) = Aa m (6 m ) = Aag (S m ) g + A (6 m ) (g)'' + • • • . (2.22) 

The size of the dominant nonsingular terms in eq. (2.21) at r = <5 ir is indicative of the size 
of Acjir- For the production of a color singlet X in the pp —> X and pp —» X +jet channels, 
the missing terms at NLO and NNLO scale as (plugging in the relevant color factors): 

qq -» X : { Acrj^ (<5ir) , Afij^ (<5jr)} ~ cr LO { — 10.7 Air In 5ir ,, 113.8 <5ir In 3 <5ir}, 

99 -> X : {Ac7j^(5ir) , Acrj^(5nt)} ~ cr LO {—24diRln5iR , 5765 ir ln 3 5nt}, 

gq -> Xq, qq ^ Xg : (Act^^ir) , Act^^ir)} « a LO { -22.7(Sir In<5 ir , 513.8 <5ir In 3 <5 ir}, 

gg Xg : {Ac7j^(5ir) , Acx^^ir)} « cr LO {-36 5iRln5iR, 1296 (5 T r ln 3 <5 IR }. 

(2.23) 


To estimate the impact of these terms relative to the full NLO and NNLO contributions, 
we write the full result for the cross section as 

+ ,(1) + ff (2)((M 2 + .... ( 2 . 2 4) 

47T \4ttJ 

We assume that the iv-factors at each order of perturbation theory for qq —> X and 
qq Xg, qg —> Xq processes are 10%, so a^/a LO ~ 10 n . For gg —> X and gg —> Xg 
processes, we assume the JL-factors are 30%, so that /a LO « 30 n for these cases. These 
factors roughly scale like the prefactors in eq. (2.23). Hence, a rough estimate of the relative 
size of the missing terms at each order is given by 


(jilt) 

(jt 1 ) 


a <5ir In <5ir , 


(frit) 

<t(2) 


a Or In 3 Or . 


(2.25) 


The dependence of these corrections on Or is plotted in figure 1, where we take a between 
1/3 and 3. The dashed line shows the known exact NLO result for thrust. This implies 
that when working to NNLO, we need Or < 10 -3 — 10 -4 to have a reasonable < 0(10%) 
determination of the ct 2 NNLO contribution to the cross section. For typical applications 
with Q ~ 0(100 GeV) this implies that Ts < 0.1 — 0.01 GeV. To the extent that the NNLO 

'In principle, subleading logarithmic terms can also be numerically important due to large numerical 
prefactors, especially for moderate Ts values. However, for small enough Ts values, the leading logarithmic 
terms are a sufficient estimate. 
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Figure 1. Estimated size of the missing nonsingular terms below r = <5 ir as a fraction of the full 
correction at NLO (blue band) and NNLO (orange band), see eq. (2.25). The dashed line shows 
the known exact result for thrust. 

terms are only a small part of the total cross section (as is the case for Drell-Yan, for 
example), a larger error on the NNLO terms might be tolerable. However, we stress that 
these estimates can only serve as an indication, and in practice one should carefully test 
the size of missing corrections, for example by studying the Jir, dependence as discussed in 
section 4.3. 

An important comment concerns the fact that it is in principle possible and straight¬ 
forward (though perhaps tedious in practice) to derive subleading factorization theorems 
for IV-jettiness and other observables using SCET. These can then be used to systemati¬ 
cally determine the next-to-singular 0(t) corrections and include them in the same way 
in the subtractions. This would substantially reduce the size of the missing nonsingular 
corrections by one power of 5 ir. A complete factorization theorem at subleading order for 
a single-jet process has been derived for semileptonic heavy quark decays in ref. [75]. For 
recent work in this direction for thrust in e + e _ see e.g. refs. [76-78]. 

A second important aspect concerns the required numerical precision in a practical 
implementation. For both 7)v-slicing and differential 7jv-subtractions, the full QCD and 
singular cross sections are probed in regions of phase space with Tn > Ts, where there 
are significant numerical enhancements due to the nearby IR singularity at Tn = 0. For 
<5ir ~ 10 , the cancellations between the full QCD and singular Tn distributions can 

easily reach the 0(1O 4 ) level and only increase as <5 ir is lowered further. Getting a result at 
O(10~ k ) relative numerical precision in this case demands at least an O(10~^ +4 ^) relative 
numerical precision in the evaluation of the squared QCD amplitudes. 

For 7Ar-slicing, the numerical cancellations only happen after the Tn integration, which 
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means that in the worst case the Tn integral itself may have to be carried out to the same 
high precision. In practice, this will strongly depend on the process and the chosen Ts, 
since the numerical cancellations actually happen between the two terms in eq. (2.17) rather 
than the the last two terms in eq. (2.18). In any case, using Monte-Carlo integration to 
determine the integral of the unsubtracted full result down to Tn > T very accurately 
requires very high statistics and good phase-space sampling. Since NLO codes are usually 
not designed for this purpose, this strongly limits how low Ts can be taken. 

For the 7]v-subtractions, the QCD amplitudes in the integrand still require the same 
high numerical precision at small Tn to obtain an accurate result for the nonsingular spec¬ 
trum. However, since the cancellations now happen already at the integrand level, the 
Tn integration itself has to be carried out only to the nominal C^ICD^) relative precision. 
Hence, the statistical requirements on the Monte-Carlo integration of the nonsingular spec¬ 
trum in eq. (2.19) are much more modest compared to the 7rv-slicing. This also means 
that Ts can now be taken as low as the numerical precision in the integrand allows. The 
main nontrivial requirement now is that one must be able to sample phase-space for fixed 
Tn, which we discuss further in section 4. 

3 iV-jettiness Subtractions 

In this section, we now specify Tn to be IV-jettiness and explicitly construct the IV-jettiness 
subtractions. We first discuss the Born kinematics and the definition of IV-jettiness in 
section 3.1. In section 3.2 we review the factorization theorem for the singular contributions 
in Tn and how the virtual QCD amplitudes enter into it. Then in section 3.3 we explicitly 
write out the Tn subtractions at NLO and NNLO. Finally, in section 3.4 we discuss how 
the subtractions can be made more differential and thereby more local. 

3.1 Definition of IV-jettiness 

3.1.1 Born kinematics 

We always use the indices a and b to label the initial states, and 1,..., IV to label the final 
states. Unless otherwise specified, a generic index i always runs over a, b, 1,..., N. We 
denote the momenta of the QCD partons in the <1 >tv Born phase space by {q a , q b ',qi, ■ ■ ■, 
and the parton types (including their spin/helicity if needed) by {k u , /«*,; ki, ..., /vat}. Thus, 
&N corresponds to 

$N = { (Qa, Ha), ( qb, H b )\ (<? 1 , Ki), . . . , {q N , K N ); $£(<?)} , (3.1) 

where &l(q) denotes the phase space for any additional nonhadronic particles in the final 
state, whose total momentum is q. (For ep or ee collisions, one or both of the incoming 
momenta are considered part of (<?)•) We will mostly suppress the nonhadronic final 
state. For us, it is only relevant because it contributes to momentum conservation in <3? tv, 
which reads 

9£ + < = 9i+-" + 9& + ?' 1 - ( 3 - 2 ) 
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When there is no ambiguity, we will associate K t = i (e.g., we use f a = f Ka ), and we use 
the collective label k to denote the whole partonic channel, i.e., 

k = {«a, Kb , Ki, • ■ •, «tv} = {a, 6; 1,..., IV} . (3.3) 

We write the massless Born momenta qi as 

Qi = Ei raf, raf = (1, Hi) , |n*| = 1. (3.4) 

In particular, for the incoming momenta we have 

E a ,b = x a , b - 7p, < = (M), ^6 =(!,-«), (3-5) 

where FI cm is the total (hadronic) center-of-mass energy and z points along the beam axis. 
The x a: b are the light-cone momentum fractions of the incoming partons, and momentum 
conservation implies 

x a E cm _ = n b ■ (q\-\ -b qN + q ), x b E cm = n a • (qi H-b qN + q) ■ (3.6) 

The total invariant mass-squared Q 2 and rapidity Y of the Born phase space are 

Q 2 = x a x b E 2 , Y = \ In — , x a E cra = Q e Y , x b E cm = Q e~ y . (3.7) 

2 x b 

The complete dd>Ar phase-space measure corresponds to 

[ d ® N = YET f— — I d ®N{q a + Qb', Qi, ■ ■ ■ , QN , q) ^ d<h L (g) ^ , (3.8) 

J 2E cm J Xa x b J 2n “ 

where d<I>Ar(...) on the right-hand side denotes the standard Lorentz-invariant IV-particle 
phase space, the sum over k runs over all partonic channels, and s K is the appropriate 
factor to take care of symmetry, flavor and spin averaging for each partonic channel. 


3.1.2 Wjettiness 

Given an M-particle phase space point with M > N, IV-jettiness is defined as [50] 

M 

Tn($m) = min { q 'n Pk i > ( 3 - 9 ^ 

ti 8 1 Qi J 


where i runs over a, 6,1,..., N. (Here we use a dimension-one definition of Tn following 
refs. [52, 62].) For ep or ee collisions, one or both of the incoming directions are absent. 
The Qi are normalization factors, which are explained below. The p^ are the M final-state 
parton momenta (so excluding the nonhadronic final state) of The qi in eq. (3.9) 

are massless Born “reference momenta”, and the corresponding directions n* = qi/\qi\ are 
referred to as the IV-jettiness axes. For later convenience we also define the normalized 
vectors 



(3.10) 
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The qi are obtained by projecting a given onto a corresponding Born point 
For this purpose, any IR safe phase-space projection can be used. That is, 
in any IR singular limit where &m —•► the Born projection has to satisfy 

< hjv( < f > .M t $tv) t ; (3-11) 


including the proper flavor assignments. In particular, for M = N, we simply have 
&n{&n) = &N and so qi = pi, which implies Tn(®n) = 0. For M > N + 1, there is 
always at least one pk that cannot be exactly aligned with any of the qi, which means 
that Tn(&m) > 0. The minimization condition in eq. (3.9) ensures that for each pk the 
smallest distance to one of the qi enters the sum, which together with eq. (3.11) implies 
that Tn(&m —> < l>Ar) —> 0. Hence, IV-jettiness satisfies all the criteria of an IR-safe IV-jet 
resolution variable given in eq. (2.4). 

Some examples of suitable Born projections are discussed in section 3.1.3 below. Al¬ 
though the precise procedure to define the Born projection and the qi is part of the definition 
of IV-jettiness, it is important that it does not actually affect the singular structure of the 
7)v-differential cross section. Different choices only differ by power-suppressed effects, as 
explained in ref. [50] , which means the precise choice only affects the nonsingular contribu¬ 
tions. Hence, constructing the singular contributions and the subtraction terms does not 
actually require one to specify the Born projection, as they are constructed in the singular 
limit starting from a given 4 > 7 y 8 This fact provides considerable freedom in the practical 
implementation, which we will come back to in section 4. 

The singular structure of Tn is determined by the minimization condition in eq. (3.9) 
and the choice of the Qi. The minimization effectively divides the 4phase space into N 
jet regions and up to 2 beam regions, where each parton in is associated (“clustered”) 
with the qi it is closest to, where the Qi determine the relative distance measure between 
the different qi. We can then rewrite eq. (3.9) as follows, 


Tn = Y j V n 

i 


with 


M 


k=i 


2g ? ; • p k 

Qi 


n»( 


Qj • Pk 
Qj 


Qi'Pk \ 

Qi ) 


(3.12) 


where the T'fa are the contributions to Tn from the ith region. 

The Qi can be chosen depending on the Born kinematics in 4 >jv (subject to the con¬ 
straint that the resulting distance measure remains IR safe). A variety of possible choices 
are discussed in detail in refs. [52, 62], An “invariant-mass” measure is obtained by choos¬ 
ing common Qi = Q. In this case, the sum of the invariant masses of all emissions in 
each region will be minimized. A class of “geometric measures” is obtained by choosing Qi 
proportional to E,j, which makes the value of Tn itself independent of the Ei, i.e., 

Qi — 2 PiEi Qi — pi * — Pi^i Pk 5 (3.13) 

^ i 

8 In this regard, the T/v-subtractions are FKS-like, namely they are intrinsically a function of the Born 
phase space 4>jv and an emission variable, which for us is Tn, as opposed to starting from a given 4>>jv+i 
point. 
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where the pi are dimensionless numbers which determine the relative size of the different 
regions. In this case, the sum of the small light-cone momenta, ni ■ Pk, of all emissions 
relative to their associated IV-jettiness axis are minimized. 

The singular structure of the cross section does explicitly depend on the distance 
measure. When discussing the singular contributions in the next section, we will keep the 
Qi arbitrary, thus enabling various choices to be explored using our results. [As discussed in 
ref. [50], one can generalize IV-jettiness further to use any IR-safe distance measure di{p k ) 
in eq. (3.9), which has been used for example in the application to jet substructure [79, 80]. 
For our purposes, the canonical form di(p k ) = qi ■ Pk is suited well, because the simple 
linear dependence on pk simplifies the theoretical analysis and computations.] 

3.1.3 Example Born projections 

To construct a generic Born projection, it suffices to use any IR-safe jet algorithm to cluster 
the M- parton final state into N jets with momenta P % . One can then define massless final- 
state = Eirf by taking (i = 1,..., N) 

P w 

Hi = -p- with Ei = Pf or E) = |P,;| or 2E* = P{ + \P%\ ■ (3-14) 

\Pi\ 

where any of the choices for Ej can be used. To ensure that the total transverse momentum 
in the Born final state adds up to zero, one can then for example boost the hadronic 
system or recoil the leptonic final state in the transverse direction. Finally, the initial- 
state momenta q a and (ft, which always lie along the beam directions as in eq. (3.5), are 
determined by momentum conservation from eq. (3.6). 

When using a geometric measure as in eq. (3.13), the canonical way to determine the 
IV-jettiness axes n* is by an overall minimization of the total value of 7 tv- Up to NNLO 
the relevant cases are M = N + 1 and M = N + 2, i.e., one and two extra emissions, in 
which case the overall minimization to find the IV-jettiness axes is still fairly easy to work 
out explicitly. 

Let us take Pi = 1 for simplicity and consider the case of hadron-hadron collisions, 
such that we have N jet axes plus the two fixed beam axes n a = ±£. When M = IV +1, it 
is easy to see that IV — 1 axes must be aligned with IV — 1 of the pk momenta. For the last 
axis, there are two possibilities, and the one which gives a smaller 7 tv is selected: Either it 
is aligned with one of the two remaining pk (this occurs if the last p k momentum lies close 
enough to one of the beam directions), or it lies along the direction of the sum of the two 
remaining pk- The appropriate expression for Tn for M = N + 1 is then: 

M 

Tn = y~](Ek - \p k \) + mini min {\pj\ - \p z A}, min { \pj | + \p k \ - \pj + p k \} } . (3.15) 

k=1 

The first term in the overall minimization corresponds to the first case above (extra emission 
clustered to the beam), whilst the second term corresponds to the second case (extra 
emission clustered to a jet). 

When M = N + 2 there are two extra emissions. Now, IV — 2 axes will always be 
aligned with IV — 2 of the p k momenta, and there are four possible cases how the remaining 
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two axes can be chosen based on the remaining four p k . The appropriate expression for 
Tn for M = N + 2 is 

M 

Tn = Ej - \pj \) + min| ™ n ^{|pj| + \p k \ - \p z \ - \p z k \}, (3.16) 

3 = 1 

. mjn {\pj\ + \p k \ + \pi\ - \pj +p k | - |pf|}, 

min I Ip,'| + |pfc| + |p| - |Pj +Pk+Pl\}, 

,, mhl + l^fel + \Pl\ + \Pm\ ~ \Pj+Pk\ - \Pl+Pm\}\ ■ 

The first term in the overall minimization corresponds to both extra particles being clus¬ 
tered to a beam direction. The second term corresponds to one particle being clustered to 
a beam, and two particles being clustered together in a jet. The third term corresponds 
to clustering three particles together in a jet, and the final term corresponds to clustering 
two sets of two particles into two separate jets. In all cases the remaining jet directions are 
set by the remaining unclustered p k momenta. 


3.2 Factorization in the singular limit 
3.2.1 Factorization theorem 

We start by writing the IV-jettiness singular cross section differential in and all indi¬ 
vidual 7 ”y contributions, 


dosing ( X ) 

d T n 

dcr sin s($7v) 

d Tn 


x(9h) 


llM 


d T n 

d U sin S($7v) 

Jd7^d7#---dT/ 


Tn~Y. T n 


(3.17) 


The factorization of the JV-jettiness cross section in the singular limit [for the linear mea¬ 
sures defined by eq. (3.9)] was derived in SCET in refs. [50-52], It takes the form 


d <r sin s($jv) 

dr^dr^---dr/ 


r N 


dta B a (t a , x a , p) / d t b B b (t b , x bl p) 


I [ / dsj Ji(si,p ) 


L i=l ' 


(3.18) 


& p) s K (r^ - ...,r/ - p^j c($ N , p) 


The first argument(s) of the beam, jet, and soft functions Bi , J*, and S K determine the 
contributions to the 7jy from the respective collinear and soft sectors. The beam function 
B a (t a ,x a ,p) contains all collinear emissions (virtual and real) from the incoming parton 
a, and depends on the parton’s flavor n a and light-cone momentum fraction x a . The jet 
function Ji(s, p) contains all collinear emissions from the outgoing parton i, and depends on 
the parton’s flavor The soft function S K contains all soft emissions between all partons 
and depends on the directions g*. It is a matrix acting in the color space of the partonic 
channel k. More precisely, it acts in the color-conserving subspace of the full color space. 
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The hard Wilson coefficient C($n,h) is a vector in the same color space, and CT($jy,/i) 
is its conjugate (see below). It contains the QCD amplitudes for the IV-parton process and 
depends on the full IV-parton phase space dw. 

All functions in the factorized cross sections have an explicit // dependence (due to their 
nonzero anomalous dimensions). This g dependence exactly cancels between the different 
functions at each order. The remaining internal //-dependence is the usual one due to the 
running of a s (/x) which cancels up to the order one is working at. In the general case, the 
// dependence is used to resum the logarithms of Tn to all orders in a s at a given order 
in logarithmic counting. For our purposes, we require the strict fixed-order expansion in 
a 3 (/i) at NLO and NNLO. 

We note in passing that starting at N 4 LO, the partonic QCD cross section receives 
a contribution from noncancelling Glauber modes in graphs with the same structure as 
figure 5 in ref. [81]. Such contributions are not reproduced by eq. (3.18). However this is 
far beyond NNLO, which is the level we are concerned about here. 

3.2.2 QCD amplitudes and color space 

The hard coefficients contain the virtual IV-parton amplitudes from QCD. They 

formally arise as the matching coefficients from QCD onto SCET. How this matching is 
performed in practice for generic processes using QCD helicity amplitudes is discussed 
extensively in refs. [82, 83] (see also refs. [84-87]). We refer the reader there for details 
and only summarize the features relevant for our discussion here. The important point 
is that when working in pure dimensional regularization with MS, the coefficients C(<&jy) 
are given by the infrared-finite part, Mfi n , of the full IV-parton QCD amplitude after UV 
renormalization. 9 Hence, we have 

C a “- a "($ N ) = -i-4fin "“' v ( < hiv), (3.19) 

where we have explicitly written out the color indices {a a ,..., ajv} of the external partons. 
(All remaining dependence on external helicities and momenta are contained in Tat.) 

The color indices {a*} span the full color space for the partonic channel n. We can now 
pick a complete basis of color structures T^ a ' aN , which span the color-conserving subspace. 
(For practical purposes, the basis can be overcomplete and does not have to be orthogonal.) 
For example, for k = gqq the color-conserving subspace is still one-dimensional, since the 
only allowed color structure is T aa13 = (T“^). For k = ggqq, one choice would be 

f aba-g = ^ T a T b )a _, , tr[T“T 6 ] S a p ). (3.20) 

Given a basis 7T“'" ajv , we write the hard coefficients in this basis as 

C a °- a N($ N ) = J2T“ a - aN C k {<f> N ) = f a “~ aN (3-21) 

k 

u The UV renormalization scheme must be the same for all functions appearing in the factorized cross 
section. The explicit results we give all use conventional dimensional regularization (CDR), which requires 
the QCD amplitudes to be renormalized in the CDR or’t Hooft-Veltman (HV) scheme. 
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This is in one-to-one correspondence to choosing a particular color decomposition for the 
IV-parton amplitude, and so the coefficients C are directly given by the IR-finite parts of 
the color-ordered (or color-stripped) amplitudes. The precise form of the amplitude’s color 
decomposition is irrelevant for our discussion and any convenient color basis can be used. 
The conjugate O' of the vector C is defined by 

C ] = E C* 0La '" 0 ‘ N T aa '"°‘ N = (C*) t T k , (3.22) 

Ota-'OIN 

where the superscript T denotes the transpose and 

T K = (fa a -a N yfa a -a N , ( 3 . 23 ) 

Ota — OtN 

is the matrix of color sums for the basis chosen for the partonic channel k. The typically 
used color bases are not orthonormal, in which case T K is not equal to the identity operator 
1 K and C' is not just the naive complex conjugate transpose of C. We then have 

|(?(<M 2 = C\$ N ) C{$ N ) = J2 |An(^)| 2 • (3.24) 

color 

3.2.3 Leading order 

It is instructive to see how the LO cross section arises from eq. (3.18). At LO, we have 

Ji°\ s i M) = S(s), 

= 5(t) f a (x,(J, F ) , 

S^\k a ,... ,k N ,{sij },n) = 1« Jj5(fei) , (3.25) 

i 

where the LO soft function is the identity operator in color space, 

1 K = 5 aapa ■ ■ ■ 5 aN ? N . (3.26) 

Plugging this back into eq. (3.18) we get 

N = fa fb C m ($N) 1^ (0) (<M 1] 8{U) 

= fa fb E |^(^)| 2 ]n^)-^)II^)- ( 3 ‘ 27 ) 

colors i i 

Equation (3.17) then reproduces the LO cross section as in eqs. (2.5) and (2.12). 


3.3 Single-differential subtractions 

We now project onto the single-differential IV-jettiness Tn- Equations (3.17) and (3.18) 
yield 


da sin z($ N ) 


d T n 


r N 

— I dt a B a (t a , x a , /r) / d tfo Bfrftbi Xbi 11 / dsj i/j(sj,//) 

J J L i=1 J 


(3.28) 


x C\d> N , fi )S K (T N -h--±-'£g, ) C(<S> N , ,i) , 

i =1 


N 
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where the single-differential soft function is the projection of the multi-differential one 
appearing in eq. (3.18), see eq. (A.21). We expand this singular contribution to the N- 
jettiness cross section as [cf. eq. (2.10)] 


d£T sin S($Jv) 
d T n 


C- 1 ($ n ,Z)6(Tn) + ^C n ^N,0 


n> 0 


(3.29) 


2 ) 71—1 




m> 0 L 


n =0 


«s(^) 

47T 


The C n (r) are the usual plus distributions defined in eq. (2.8). Here we explicitly denote 
the dependence of the subtraction coefficients £, n) on the renormalization scale 

/j. The individual coefficients also depend on the arbitrary dimension-one parameter £, 
which drops out exactly in the sum of all coefficients at each order in a s . (In section 2.2 we 
used £ = Q.) Finally, the coefficients also depend on the IV-jettiness measures Qi, which 
we suppress for simplicity. 

To determine the subtraction coefficients, we simply expand all the functions in the 
factorization theorem eq. (3.28) in terms of a s (/i), 

m> 1 

= S(t ) f a (x,HF ) + , 

rri> 1 

= i K 5(k) + Y ^(MsYU^^y 71 , 

rri> 1 

= <? (0) (<W) + £ cM (<W )(^)) m , (3.30) 

m> 1 

plug these back, and collect all contributions to each order in a s and each power in lnT/v- 
Explicit results for the jet, beam, and soft functions through O(o%) are given in Ap¬ 
pendix A. 


3.3.1 NLO subtractions 

At NLO, the differential subtractions require the subtraction coefficients Cq 1 ^ and C^\ 
which are the coefficients of the 1/Tn and (InTivO/T/v contributions. They are given by 
(with n = 0,1), 


cW(^,£, M ) = |c(°)(<f> Jv ^)| 2 


fa(x a , /If) fb(x b , Hf) 


E 4 ’(— 

i=l 




fi B a,n (x a , Ab HF, fb(Xb, VF) + fa{x a , Hf) (x b , fj,, (Ip, ^)| 

+ fa(x a , fxp) f b (x b , lip) Ct(°)(^, /i) Sjti ({ft}, i) C (0) (<fw, m) • (3.31) 
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The jet function contributions in the first line effectively correspond to collinear final- 
state subtractions, while the beam function contributions in the second line effectively 
correspond to collinear initial-state subtractions. The soft function contribution in the last 
line effectively corresponds to a soft subtraction. 

As explained in section 2, the coefficient C_i determines the integrated subtractions 
plus the virtual contributions, which becomes obvious when choosing £ = T Q s- At NLO, 
we have 




fa{x a , ti F ) f b (x b , n F ) (Ct(Dc(0) + cKo)^!)) ($JVj ^ 

r N op'- 

+ |C ( 0 ) (d>AT,(r )| 2 f a (x a , Hf) fb(x b , fJ, F ) 

L ?;=i ^ ' 

+ B a~ 1 (xa, /b VF, fb(x b , Hf) + fa{x a , Hf)B^]_ ± (x b , (J., /Ip, 
+ fa(x a , Hf) fb(x b , Hf) C t(0) ($N, I*) S^l _! ({&}, C (0) ($jv, fJ,) . 


Qb^\ 
fi 2 J 

(3.32) 


The hrst line contains the IR-finite virtual one-loop amplitudes in The remain¬ 

ing lines effectively correspond to the integrated collinear and soft subtractions. The NLO 
beam, jet, and soft function coefficients, B^n(x, /x, fi F , A), ,J- ] J (A), and si ] ,n({q t }. A) ap¬ 
pearing in eqs. (3.31) and (3.32) are all known and are collected in Appendix A. The PDF 
factorization scale dependence only enters via the beam functions and the PDFs. 

One can see that the structure of the subtraction terms has a close resemblance with 
FKS subtractions. The important difference is that here one does not have to divide up 
phase space in order to individually isolate all possible IR singular regions. Instead, all 
the singular regions are projected onto the single variable Tn- An analogous phase-space 
division for the soft emissions now happens in the calculation of the N-jettiness soft func¬ 
tion. It is also important to note that there are no overlaps (i.e. double counting) between 
the soft and collinear subtraction terms. In principle, such overlaps can exist and must be 
removed, which in SCET corresponds to removing so-called zero-bin contributions [88]. A 
nice feature of IV-jettiness is that all such overlap contributions automatically vanish in 
pure dimensional regularization at all orders in perturbation theory. 


3.3.2 NNLO subtractions 

For simplicity of the presentation, we define the abbreviations 


tM — tM / Qj£ \ 
°i,n ~ u i,n V (X 2 / 



oM _ dM 
D z,n ' C 2,n 


cM\ 
l^F•> ^2 J ? 


& m) =£(m) 


fj — fi(Xi , [J-f) ; 

(3.33) 


where we use roman letters (B, J, S, C, f) to avoid any confusion with some of the coefficients 

( 2 ) ( 2 ) 

listed in Appendix A. The NNLO coefficients J) and B) ’ as well as the soft function 
W2) 

coefficients are all known analytically, see Appendix A. 
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The IV-jettiness soft function describes how the soft radiation is split into the different 
IV-jettiness regions. Obtaining the two-loop soft constant S y K _ 1 (the coefficient of the 8(k), 
see eq. (A.23)) is the remaining principal challenge. It is known analytically for processes 
with two external partons, see eqs. (A.34) and (A.35), but it is currently unknown for 
generic IV-jet processes. It can however be determined numerically by extending the NLO 
calculation in ref. [52] using existing NNLO results. A procedure to do so has been outlined 
recently in ref. [89], where numerical results for 1-jettiness in pp collisions were presented. 

We conveniently denote the genuine m-loop contributions from jet, beam, and soft 
functions to the Cn^ as 


x(m) = |c(°)| 2 (i a f 6 ^ j£? + Bg f 6 + f a + f a f 6 Ct(°) sW C<°>. (3.34) 

' i =1 ' 

Using this notation, we write the two-loop cross terms related to real-virtual contributions 
and involving the one-loop virtual amplitudes in CL 1 ) as 


N 


4 1+1) = {C m C^ + Ct^cW) f„ f 6 £ Jg + Bj) f 6 + f Q B« 




2=1 


+ f a f fe (C t(0) c (1) + cWsWcW) 


(3.35) 


Finally, the cross terms with two one-loop coefficients of jet, beam, or soft functions from 
the associated C n <8> C m convolution are denoted as 


N 


N 


N 


4S 1 ’ = |c (0) l 2 (fui, y . J 2 J 7L + + B S B S 

^ i<j = 1 2=1 2=1 

N 

+ fa fb J S C t(0) S&> C ( °) + B« f b 3t(°) SW C<°> + f a Bg s« C (0) . 

2=1 

(3.36) 

With these definitions, the NNLO subtraction coefficients read 

ci 2) (<f> N ,z,fi) = xi 2) + x £ +1) , 

4 2 ) ($ w , c , m ) = 4 ° + k4T ] + X IT ] ) - 


C< 2) ($*,£,/*) = x (2) + a{ 1+1) + 2 A^ 0 +1) - , 


cS 2) ($iv,^)=4 2) +4 i+i) - y(4; 1 +i )+<o +1 ))+2c 3 < i) , 

= f„f 6 (c t(0) c (2) + c t(1) c (1) + c t(2) c (0) ) 

+ x® + x B+,) - L x a+i)+cafx'y 11 +x B+I) ) - . 


(3.37) 

(3.38) 

(3.39) 

(3.40) 

(3.41) 


(2) 

The S(Tn) coefficient C_{ again corresponds to the integrated NNLO subtraction piece and 
contains the full IR-finite 0(a 2 ) virtual IV-parton amplitudes in CL 2 ). 
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k 

v k 00 

y k 01 = v k w 


-1 

— 7T 2 /6 

C 3 

O 

CO 

CO 

1 

0 

0 

— 7T 2 /6 

2Ca 

1 

2 

0 

— 7T 2 /3 

2 

0 

3/2 

0 

3 

0 

0 

1 


Table 1 . Coefficients V/™ 71 for the convolution C m <g> £„ according to eq. (3.42). 

The constants multiplying the Xn,m in equations (3.37)-(3.41) are the coefficients 
V™ n arising in the convolution C m . <g> C n , 


/ m+n+l 

d T'C rn (T-T')C n (T , ) = VLTS(r)+ £ v r^k(r). (3.42) 

k =0 

They are given in table 1 for m,n < 1. Their expression for general m,n can be found in 
Appendix B of ref. [74]. 


3.3.3 Toward N 3 LO subtractions 

Using the notation introduced in the previous subsection it is straightforward to also write 
down the N 3 LO fV-jettiness subtraction terms. Besides the genuine three-loop terms 
according to eq. (3.34), we now have “two-loop times one-loop” cross terms Xn +2> = Xn +l ^ 
and Xn,m ^ = Xn^ as well as the “(one-loop) 3 ” cross terms Xn +l+1 \ Xn£^ +l \ and 
X^\ where the latter is associated with the convolution C n <g> C m <8> £/. 

The N 3 LO subtraction coefficients then schematically take the form 

Ccj 3 \<&]\r,€,/x) = X^ + cross terms, 


: (3.43) 

h) = Xj 3) + cross terms, 

C {3 l($ N , £, fi) = X {3 1 + f a f b + C t(3 )C ( °)) + cross terms . 


The cross terms in eq. (3.43) are a linear combination of the above listed X’s, whose 
numerical coefficients can be easily worked out by evaluating the relevant convolutions 
among the £ n <3 distributions in analogy to the NNLO case. 

For processes with only two colored external partons, so e + e _ — >qq, DIS, or pp —> color 

(3) 

singlet, analytic expressions for all X^ 0 are in fact available. This is because the three- 
loop anomalous dimensions of jet and beam functions, the PDFs, and the hard function 
are known [90-97], which also fixes the three-loop soft anomalous dimension. This means 
the complete set of 0(a 3 ) logarithmic (C n ) terms of the renormalized jet, beam, and soft 
functions are determined by their RGE. The only coefficient that is not fully known is C_{, 
associated with the integrated N 3 LO subtractions plus virtual corrections. The C/ 3 ' are 
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known from the IR-finite parts of the three-loop quark and gluon form factors [98, 99]. 
As the cross terms only involve known lower-order contributions, the only missing piece 
in eq. (3.43) then is X_-[, for which one has to compute the three-loop /r-independent 
constants of the jet, beam, and soft functions. 


3.4 Constructing more-differential subtractions 


As mentioned already, the 7j\r-subtractions we have defined thus far are nonlocal, in the 
sense that all the singular regions are projected onto the single variable Tn and the subtrac¬ 
tion acts only after the corresponding phase-space integrations. It is conceivable that in 
order to improve the numerical stability or convergence of the NNLO calculation one might 
wish to use a more local subtraction - indeed, many of the available NNLO subtraction 
schemes utilize highly local subtraction terms. 

In our approach, it is straightforward, at least conceptually, to progressively increase 
the locality of the subtractions. All one needs to do is split Tn up into further IR-safe 
observables that cover the phase space and which are sensitive to emissions in different 
regions, and/or introduce further observables that resolve the nature of emissions, e.g. 
allowing one to discriminate between double-real and single-real(-l-virtual) emissions in a 
given region. The subtraction is then given by the singular cross section differential in all 
of these observables. In practice, this requires the relevant factorization theorem for this 
more-differential cross section. 

Let us demonstrate how this works for a simple example. The factorization theorem 
in eq. (3.18) is already differential in the individual A r -jettiness contributions T/y. For 
simplicity, we take IV = 0 and consider the X + 0 j NNLO cross section. In this case, 0- 
jettiness (aka beam thrust) effectively splits the event into two hemispheres (beam regions) 
a and b, whose A r -jettiness axes are defined by the beam directions. The total 0-jettiness is 
given by 7o = 7/° + Tq, where T { °' and Tq are the contributions from the two hemispheres 
[cf. eq. (3.12) and its discussion]. 

Following the procedure in section 3.3 we can use the total To to construct a subtrac¬ 
tion. However, instead of taking the sum, we can also consider T a = Tq and % = 7// 
separately, and perform the subtraction differential in both of these observables. Each 
of them is then sensitive to a subset of the singular regions, namely, collinear (and soft) 
emissions closer to beam a will only affect T a , whilst emissions closer to beam b will only 
affect Tb- 

Following the logic of section 2.3, we first write down the appropriate formula for the 
corresponding double-differential phase-space slicing: 


rTs r 

a(X) = d T a 

Jo Jo 


/ o d7ad7fe J j- s Jo dT a dTb 


rTs 


+ I d Ta [ d T b -^-t + / d Ta f d Tb + 0(5 m ) . (3.44) 

Jo Jt s 0.1 a d lb Jjg Jt & Ol a Olb 


Here, we substitute in the double-differential singular cross section when both T a and Tb 
are below the IR cutoff Ts, which is correct up to 0(<5ir). Having either T a or % nonzero 
requires at least one additional emission, so the remaining three regions only require an 


- 23 - 






Figure 2. Division of the T a ,Tb phase space for double-differential TV-jettiness subtractions. 


NLO calculation. Of course, there are singularities in the second term as T b —> 0 with 
nonzero T a (and similar singularities in the third term when T a —> 0 with nonzero 7ft), but 
these are handled as part of the NLO calculation. 

Performing the slicing method using both T a and 7ft has no clear advantage over the 
slicing method using To alone, as in both methods one basically removes a small region 
of size Ts around 7o = 0, and handles it using the singular cross section. However, let us 
rewrite eq. (3.44) as a subtraction by adding and subtracting the singular cross section for 
the shaded region in figure 2, arranged in the following way: 


ct(X) = cr sing (X, T a < T oS ,T b < r off ) 


rTos 

+ / d T h 

JTs 
rT 0 s 

+ / d Ta 

Jt s 

rT 0 f[ rT 0 ff 

+ d Ta dTft 

Jt 5 JTs 


MX, Ta < Ts) da si ^(X, Ta < Ts) 
dTft dTft 

MX, Tb < Ts) da si ^(X, % < Ts) 
d Ta d Ta 

d a{X) dcr sing (X) 


d T a dT b dT a dT b 


+ d Ta dTb 


da(X) 

dTa dTb 


[1 - e(T a < Ts) 9(T b < T 0 ff)] + 0(Sm ). (3.45) 


This equation is the two-variable analogue of eq. (2.19). The parameter Toff controls again 
where we turn off the subtraction, and the dependence on it precisely cancels between all 
contributions. The total cumulant in the first term contains the two-loop virtual correc¬ 
tions together with the corresponding integrated subtraction terms. The cross sections in 
the second and third terms are differential in one of the variables and integrated in the 
other. Since one of the variables is nonzero, while the other is integrated, they require an 
NLO calculation with one additional resolved emission. These terms contain all the real- 
virtual contributions and the singular cross section acts as the corresponding real-virtual 





















subtraction. The fourth term involves the double-differential cross sections and since both 
variables are nonzero only requires a LO calculation with two resolved emissions, one in 
each hemisphere. The double-differential singular cross section then acts as the correspond¬ 
ing double-real subtraction, which is point-by-point in both T a and Ti,■ (Contributions with 
two real emissions in the same hemisphere are part of the NLO calculations in the second 
and third terms.) Hence, by considering separately T a and %, one is able to disentangle 
different real-virtual and double-real contributions and also make the subtractions more 
local. The price one has to pay is that the double-differential singular cross section in the 
last term requires the double-differential NNLO soft function, which is more complicated. 
(For the beam and jet functions this requires no additional effort.) 

A further important point to make is that T a and Tb are defined such that requiring 
T a > Ts or Tb > Ts forces the corresponding emission to be in hemisphere a or b. At NLO, 
there is only one real emission, so only one out of T a and Tb can be nonzero. Then, the 
double-differential subtraction essentially splits the 7o-subtraction into two pieces, acting 
in the two hemispheres. At NNLO, this splits the real-virtual contributions into the two 
pieces in the second and third lines of eq. (3.45). If this is undesired, one can instead 
consider the two variables T m in = min{71,,, T,} and 7jnax = max.{T a ,%}- This effectively 
folds the phase space in figure 2 in half along the diagonal where T a = T, and combines 
the second and third terms in eq. (3.45) into one. 

Now let us return to the general case with N partons in the Born process. Then 
there are N + 2 contributions T^,T^,T^, ■ .. ,T^, and one can consider the subtraction 
separately in all of them. At NLO, only one of them can be nonzero, while at NNLO at most 
two of them can be nonzero. This means that there will be many different contributions, 
where in each contribution only one or two of the T^ are differential and nonzero, while 
all the others are integrated over. Each T^ can only be nonzero when the corresponding 
emission is in the ith IV-jettiness region. Hence, if desired, using the individual 7~{- as the 
resolution variables automatically yields a division of phase space into different singular 
regions around each of the N partons, very similar to the phase-space divisions encountered 
in traditional local subtraction methods. On the other hand, if the proliferation of phase- 
space regions is undesired, one can still have the same gain at NNLO as in eq. (3.45) by 
considering two combinations of all 7jy, e.g. the minimum and maximum nonzero TJy, or 
the sum of all together with the sum of all but the largest T^- 

Instead of or in addition to splitting Tn into its different components, one can also in¬ 
crease the locality of the subtraction by performing it differentially in both Tn and another 
independent IV-jet resolution variable. For example, one could look into each IV-jettiness 
region i and compute the scalar sum of transverse momenta with respect to the corre¬ 
sponding IV-jettiness axis Exi = ^2k&\PTk\, performing the subtraction also differential 
in the Ex%■ Doing so resolves part of the radiation phase space, which would otherwise 
be integrated over when considering only Tn by itself. For the X + 0 j case one could 
for example consider 7o together with the transverse momentum px of the color-singlet 
final state X. The relevant factorization formulae differential in 7o and px have been dis¬ 
cussed and written down in refs. [63, 100] (see also refs. [101, 102]), and the corresponding 
double-differential two-loop quark beam functions have been computed in ref. [103]. 
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We have discussed several options how to extend the single-differential IV-jettiness sub¬ 
tractions, but of course this is not an exhaustive list. Constructing such more-differential 
subtractions requires the appropriate singular cross section differential in all of the chosen 
jet resolution variables, and in order to experience the maximum advantage in terms of con¬ 
vergence, these differential cross sections should reproduce the correct singular behaviour in 
all of the relevant singular kinematic regimes. The factorization of multi-differential cross 
sections in SCET accurate in all relevant kinematic regimes is a topic that has received 
much interest recently, see e.g. refs. [52, 61-65], and it would be interesting to apply this 
work to the issue of calculating NNLO QCD cross sections. 


4 Practical Considerations and Implementation 

In this section, we discuss in more detail how the singular cross section in eq. (3.29) 
can be implemented in practice as a subtraction term following our general discussion in 
section 2.3. We first discuss the NLO case in section 4.1, where we also highlight the 
similarities to FKS subtractions, and then the NNLO case in section 4.2. In section 4.3, we 
discuss numerical aspects using the NNLO rapidity spectrum in Drell-Yan and gluon-fusion 
Higgs production as an example. 


4.1 NLO 

4.1.1 FKS subtractions 


In the notation of eq. (2.1), the cross section at NLO is given by 
(T nlo (X ) = f d ^ N B N ^ N )X(^ N ) 


+ 


d4>jv Vn(®n) Y(4>jv) + Jd<f> N+1 B N+1 (<S> N+1 )X(<Z> N+1 ) 


(4.1) 


e-s>0 


where Vjy is the IV-parton virtual one-loop contribution and Bn +i is the N + 1-parton 
real-emission contribution, 

VN(*N) = fa h + •4 I 1a4'Liv] (®A0 • 


color 


B N+1 ($ N+1 ) = f a f b ^ (^TV+l) | 2 ■ 


(4.2) 


color 


The additional a s in -B/v+i compared to Bjy is contained in ^ 4 !°) (<p 7 V+ 1 ) • As indicated in 
eq. (4.1), the limit e —> 0 can only be taken in the sum of Vm and integral over Hjv+i- 
When implementing eq. (4.1) using FKS subtractions [1-3, 104, 105], the cross section 
is obtained as follows: 


a NLO (X) = J d$ N |(Bjv + Vg)($ N )X($ N ) 

+ E/ d<J> rad UB k N+1 X)(^ k N+1 )-S k N+1 (^ N ,^ d )X^ N ) \+0(6 m ), 


'< 5 ir 


vg($ N ) = 


Vn(®n) + ^2 j d^rad Sn + i(®Ni ^rad) 


(4.3) 


£->0 
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Here, the phase space is first sampled over <3 ?tv- For a fixed <&n point, one then further 
samples over the radiation phase space H> rac j, where the sum over k runs over all the different 
IR-singular regions. The real-emission contribution and measurement (.Bjv+i-^X^iV-i-i) = 
-B/v+i(3 ?jv+i)-^(^jv+i) are evaluated at a constructed point 4>X+i = ^v+i(^Nj 4> ra d)- The 
superscript k on B^ +1 indicates that Bn+i is divided up between the regions in such a 
way that it is precisely reproduced in the sum over all regions. The phase-space map &n+i 
and the subtraction terms S//+i are specific to each singular region. The S^ +1 are directly 
constructed in the singular limit, meaning they are functions of 4? at and Trad only, and in 
particular do not depend on the actual map &%+i- In practice, there is again a tiny IR 
cutoff 5 ir required on the <h ra d integral due to limited numerical precision and the fact that 
B k N+1 and S^ +1 each individually diverge. The subtracted virtual, V$, contains the finite 
remainder after combining the virtual contributions with the integral of the subtractions 
and cancelling all 1/e IR poles. 


4.1.2 T/v-subtractions 

As discussed in section 2.3, the full cross section for X at Tn > 0 only requires a lower-order 
calculation. At NLO, we need its LO expression given by 

LO p 

= / d<f>AT+i (Bn+iX)($n+i) 5[Tn - Tn(®n+i)] , (4.4) 

T/v>0 J 


d a{X) 
c\T n 


where it is obvious that this is a LO quantity. 

Since the subtractions are used up to the upper cutoff Tn < T Q s, as seen in eq. (2.19), it 
is most convenient to set £ = T 0 s in the subtraction coefficients. For the singular spectrum 
at Tn > 0, we can simply drop C_i and replace C n (r) —y In n (r)/r. The subtraction terms 
at NLO are then 


<r sing (<I>Ar, Tos) = [ T ° S d Tn daS ™ S (® N) = C-^jv.Tofr), 

Jo d In 

1 


da sing ($jv) 


d Tn 


Tn>0 


Tn 


Co{*N, Toff )+C 1 (*N, Toff) ln(^) 


0(T n < To ff ), (4.5) 


where for convenience we included the 9(Tn < T>ff) in the singular spectrum. 
Using the above with eq. (2.17), the Tv-slicing at NLO becomes 


a NLO (X) = J d^ N a sin ^ N ,Ts)X^ N ) 

+ Jd$ N+1 (B N+1 X)($ N+1 ) 6 [T n ($n+ i) > Ts\ + 0(5 m ) . (4.6) 


This calculation is very easy from an implementation point of view, since it boils down 
to performing two LO phase-space integrals. As already eluded to, the main practical 
limitation are the large numerical cancellations between both terms, requiring the phase- 
space integrals to be evaluated to very high precision. 
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Using eq. (4.4), the differential Tjv-subtraction in eq. (2.19) takes the form 
a NLO (X) = I d<S> N a sin z(<S> N ,T oS )X(<5> N ) 

+ J <1Tn ^Jd&N+i (Bn+iX)($n+i) $[Tn — Tn(®n+i)] 

- J d$ 7 v daS ^ Jv) X ($ N )| + 0(<s m ). (4.7) 

For a numerical implementation, one must be able to solve the 5 function in the d<hjv + i in¬ 
tegral, which amounts to being able to sample over all of ‘hAr+i that gives a fixed Tn(&n+i)- 
One option to do so is to decompose <I> 7 v+i as 


®n+i = 4?at <8 > Tn <8> $Tad, d4>7v+i = d4>Ar dT/v dfl ra d , (4.8) 


where 4 >tv = &n($n+i) is precisely the Born projection used to define Tn($n+i), see 
section 3.1. The O ra d = Hrad^jv+i) contains the remaining information needed to fully 
specify 4 >at+i, which includes the continuous angular radiation variables as well as the 
discrete information about flavor, spin, and in which IV-jettiness region the additional 
emission goes. We can then rewrite eq. (4.7) as 


a NLO (x) = f d$JV J <T smg(^ Toff) X ($ n) 


+ / dT N 
Jt s 


^iS4 ra d (Bn+iX)(<&n®Tn® S4rad)- N - X($n) 


(4.9) 

+ 0(<5m). 


One now samples first over 4 >at and then Tn- For fixed 4 >tv and Tn, one further samples 
over fl ra d and evaluates the real-emission contribution at the 4 >at+i point reconstructed 
from all of these. Being able to reconstruct <&jv+i(4>jv, 7)v; ^rad) is equivalent to inverting 
the Born projection. Recall however, that the singular contributions dcr smg (<I>jv)/d77v are 
independent of the Born projection. Therefore, one has the freedom to specifically choose 
the Born projection to facilitate this inversion, making it easily possible. 

Note that the O ra d integral contains a discrete sum over all IV-jettiness axis/regions. 
If one were to separate Tn into its individual components T ^ as discussed in section 3.4, 
this sum would become explicit and the single subtraction term dCT smg (4>Ar)/d71v would 
effectively separate into different subtraction terms for each region. 

We also note the close similarity of eq. (4.9) with the FKS subtraction in eq. (4.3). 
Basically, Tn <S> f2 ra d now acts as <3? ra d, while the split up into singular regions is now deter¬ 
mined by the definition of Tn- The LO piece C_\ of <7 sing supplies the Born contribution 
Bn, and the NLO piece corresponds to V^ T . 


4.2 NNLO 

At NNLO, the subtraction terms are 

<j sms ($N,T 0 g) = C-i($n,%s) , 
d a sing ($jv) 


dT N 


Tn>0 


^J2 c n ($N,T c iff) ln”(^) 9(T n < Toff) 

n =0 <n 


(4.10) 
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As at NLO, we have chosen £ = T 0 s and included the 6{Tn < 7off) in the singular spectrum. 

The full T/v-differential cross section at T/v > 0 is now needed at NLO, where it is 
given by 


^ ={ [d$N+i (Bn + iX+ Vn + iX)((!>n + i) 5[Tn — Tn($n+i)] 

d ' N T n > o U 

+ / d$ N+2 (B N+2 X)($ N+2 ) 5 [Tn ~ Tn(®N+2)] > 

■J J e-s>0 

= Jd$ N+1 1 (B n+1 X + Vg +1 X)($ N+1 ) S[T n - Tn($n+ i )] 
+ E f d< ^ad \(B k N+2 X)($ k N+2 ) 5 {T n - T N ($ k N+2 )] 


- S k N+2 ($ N+1 , ^ rad ) X{$ N+1 ) 5[T n - 7^(3* at+i)] 


(4.11) 


In the second equation we wrote it in the form of an NLO calculation with FKS-like 
subtractions, analogous to eq. (4.3), where now &% +2 = 4*^+2(^W+i 5 4> rad ). 

In general, the &n +2 map used in the N + 1-jet NLO calculation will not preserve 
Tn, that is, 7}v+ii ^rad)] / Tn(&n+i)- This means we have to be careful in 
implementing the Ts cutoff, because in order for the neglected pieces to be nonsingular, 
the cutoff must be applied on the true Tn($n+ 2 )- We can do this by treating the cutoff 
analogous to the measurement X. That is, we define the NLO calculation 

d C: ( ^ = (B N+I x + V$ + 1 X)(* N+1 ) 9[Tn($n+i) - Ts] 

dtPjv+i 

+ d^rad {(Sjv+iX)(^ +2 ) 9[T N ^ k N+2 ) ~ Ts] 

k 

- S% +2 {$ N+ 1 , $ rad ) X($ N+1 ) 6[T n ($n+ i) -Ts]}, (4.12) 

which is fully-differential in <l? 7 v+i and satisfies 



d<r(X) 

dT N 


NLO 
Tn> 0 


J d4>7v+i 


du NLO (A,^) 

d^Ar+i 


(4.13) 


Using the above together with eq. (2.17), the 7Ar-slicing at NNLO is given by 

a NNLO (X) = [d<S> N a sin z(<S> N ,T 5 )X($ N ) + [ d4>7v +1 d<7NL ° (X »^) + 0(«s IR ) . ( 4 .14) 

J J Ci4>Ar + i 


This is again quite easy to implement, only requiring a LO phase-space integral for the first 
term and a NLO calculation for the second term. The practical limitation is the achievable 
numerical precision in the NLO calculation and the ^n+i integral, which strongly limits 
how low 7 ~s can be pushed. 
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From eq. (2.19), the differential 7Ar-subtraction at NNLO takes the form 


a N NLO(X) = / d<S> N a sin H$ N ,T oS )X(<S> N ) 


+ / dT N 
JT S 


NL ° - x{ltN) 


d Tn 


7jv>0 


d T n 


+ 0{5 ik ). (4.15) 


To implement this numerically, one must be able to compute the Tn spectrum dcr(X)/d7jv 
to NLO for a given Tn, which requires to solve the <5 functions in eq. (4.11). For 4>Ar_|_i 
we can use the same procedure as at NLO together with the N + 1-jet NLO cross section 
dcr NLO (A, Ts)/ d<hjv+i> giving 


a N NLO(X) = / d^TV <! * smK (<S>N, Toff) X($ N ) + / d Tn 


dflrad 


da NLO (X,T s ) 


d^Ar+i 


d ff sing ($Jv) 

d Tn 


X(Q n ) 9{T n - Ts) 


+ 0 { 8 m ). 


<J>iVlS)7iv®^rad 

(4.16) 


This can be implemented like the NLO case in eq. (4.9), with the LO Bn+i{&n+i) replaced 
by dcr NLO (A, 75)/d<J>jv+i- The subtraction in eq. (4.16) is not completely local in Tn, since 
the 4 >at+ 2 points being integrated over in dcr NLO (X, 75)/d4>jv+i will generically not have the 
correct Tn value. However, a simple phase-space map , 2 that approximately preserves 
Tn might be sufficient in practice. 

To achieve an exact point-by-point cancellation in Tn, one also has to solve the 8(Tn — 
Tn(&n+ 2 ) constraint in eq. (4.11). This requires constructing a map 4*^, 2 f° r the IV+1-jet 
NLO calculation that preserves Tn so Tn[&n+ 2 (®n+i, 4> ra d)] = Tn(&n+i) and is equivalent 
to inverting the Born projection &n(&n+ 2 ) underlying the definition of Tn(&n+ 2 )- This 
is quite a bit more challenging than at NLO. It has been achieved in ref. [67] for a slightly 
modified version of Tn- Assuming, we have a &n +2 nrap like this, we can pull the Ts cut 
out of the NLO calculation, such that 


d a NLO (X,T 5 ) 

dTjv+i 

d<7 NLO (X) 

d4>Ar + i 


du NLO (X) 


6[Tn($n+ 1 ) - Ts] 


d&N+i 

= (B n+ iX + Vg +l X)($ N+ 1 ) 


(4.17) 



(■ B N+ iX)(<S > k N+2 ) - S k N+2 ($ N+ i, 4> rad ) X($ N+1 ) 


The differential 7j\r-subtraction then becomes 


<t nnl °(x) = j d$ N |a sin g(^,r off ) x($ N ) 


+ [ d T n 
Jt s 



dq NLO (X) 

d4>jv+i 


<J , JV(S)7jv ( S)^rad 


dcJ sirlg (4>Ar) 

d7^ 



+ 0(«J m ). 


(4.18) 


The subtraction is now fully localized in Tn, and the only nonlocality is in the dH ra d 
variables. 
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Figure 3. The absolute value of the full, singular, and nonsingular contributions to the 7o spectrum 
for Drell-Yan production. The NLO 0(a s ) corrections are shown on the left, and the pure NNLO 
O(ag) corrections are on the right. 




Figure 4. The nonsingular 7o spectrum for Drell-Yan as a function of r = To/mz- The NLO 
0(a s ) corrections are shown on the left, and the pure NNLO 0(a 2 s ) corrections are on the right. 


In a practical implementation of eq. (4.16) or eq. (4.18), the subtraction terms are 
easily evaluated and the most nontrivial ingredient is in fact the NLO calculation of 
dcr NLO (X)/d4 > Ar+i, for which one can use any existing FKS-like NLO calculation or one 
can iterate the IV-jettiness subtractions and perform it using 7/v+i-subtractions. Note that 
in all cases above the X measurement is performed inside d(T NLO (X)/d4 > jv-)-i- If the < hjV _)_2 
map preserves X , so X(&^ +2 ) = A"(4>tv+i), then it can be pulled out of the N + 1-jet NLO 
calculation. 

4.3 Example: NNLO rapidity spectrum for Drell-Yan and Higgs 

To illustrate our method with a nontrivial example, we consider the rapidity distribution 
of the vector boson in Drell-Yan production, pp —> Z/'y —> , and of the Higgs boson in 

gluon fusion, gg —> H, which are known to NNLO [39, 40, 106-110]. Since the size of the 
perturbative corrections in the two cases are very different, they provide very useful and 
complementary test cases. 

In both cases, 0-jettiness 7o is the resolution variable and all of the ingredients neces¬ 
sary to implement the 7o-subtractions through NNLO are known. (We use the geometric 


- 31 - 



















Figure 5. The absolute value of the full, singular, 
for gluon-fusion Higgs production. The NLO O(o 
NNLO O(a^) corrections are on the right. 



and nonsingular contributions to the 7o spectrum 
s ) corrections are shown on the left, and the pure 



Figure 6. The nonsingular To spectrum for gluon-fusion Higgs production as a function of t = 
To/mn ■ The NLO 0(a s ) corrections are shown on the left, and the pure NNLO 0(0%) corrections 
are on the right. 

measure with pi = 1, see eq. (3.13), which makes 7o identical to beam thrust.) The re¬ 
sults are obtained for the LHC with a center-of-mass energy of 13TeV. We always use 
CT10 NNLO PDFs [111]. We choose common renormalization and factorization scales, 
put = jip = p. with p = mz for Drell-Yan production and p = mu for Higgs production. 
For the latter we use mu = 125 GeV and work in the top EFT limit. For the Z + 1-jet and 
H + 1-jet NLO calculations we use MCFM [112, 113]. 

An important validation of the IV-jettiness subtractions is to confirm that the singular 
Tn spectrum is correctly describing the Tn —> 0 singularities of the full QCD result. 
This is done by calculating the nonsingular Tn spectrum as in eq. (2.15) as the difference 
of the full QCD and singular Tn spectra. The decomposition of the To spectrum into 
singular and nonsingular components is shown in figures 3 and 5 for Drell-Yan and Higgs 
production, respectively, where we separately show the 0(a s ) (NLO) and O(a^) (pure 
NNLO) corrections, counted relative to the LO Born cross section. (We plot the magnitudes 
of the contributions on a logarithmic scale, and the dips at large To and around 7o = 1 GeV 
are due to the spectra going through 0. The small jitters in the pure NNLO nonsingular 
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are due to numerical inaccuracies.) One can clearly see the large numerical cancellations 
between the full and singular results for small 7o, where the nonsingular spectrum is several 
orders of magnitude smaller than the full and singular spectra. 

As shown in eq. (2.19), the nonsingular spectrum is precisely the quantity that one 
integrates numerically when using the differential 77v-subtractions. The fact that the non¬ 
singular only contains integrable singularities is seen in figures 3 and 5 by its smaller slope 
toward 7o —> 0. To check explicitly that the subtractions work and the nonsingular does not 
contain any 1/Tn singularities, we consider the distribution der nons /d In r = rdcr noris /dr, 
which must go to 0 in the Tn —>• 0 limit. We plot it in figure 4 for Drell-Yan and in figure 6 
for Higgs production, again separately for the NLO and pure NNLO corrections, and using 
r = To/mz and t = To/mn, respectively. One can see that d<7 nons /dlnT —> 0 for r —> 0, as 
it must. The error bars come from the statistical integration uncertainties in the full result 
obtained from MCFM. The numerical uncertainties in the singular result are negligible in 
comparison. 

To obtain results for the NNLO rapidity spectrum, we use the simple 7o-slicing in 
eq. (4.14). As explained earlier, the missing 0(<5 ir) contributions due to the Ts cutoff 
are the same irrespective of how the subtractions are implemented. At NLO, we use 
Ts = 0.03 GeV (<5ir cs 3.2 x 10 -4 for Drell-Yan and <5ir = 2.4 x 10 -4 for Higgs) and at 
NNLO we use Ts = 0.1 GeV (<5ir k 1.1 x 10 -3 for Drell-Yan and <5ir = 8 x 10 -4 for Higgs). 
These values are at the lower end of r values plotted in figures 4 and 6, and are mainly 
limited by the MCFM statistics. 

Specifically, we use MCFM to compute the NLO cross section for 7o > Ts in the 
second term in eq. (4.14), 


d$i 


dcr NLO (Y, Ts) 


(4.19) 


J (IT i 

in bins of Y for the processes pp —> Z/'y —> £ + £~ + jet and pp —» H + jet. Since there 
are no cuts on the final-state jets other than the requirement To > Ts, for small Ts the 
calculation probes deep into the singular region and care must be taken to obtain reliable 
and numerically stable results. This is combined with our own implementation of the 
NNLO singular cross section for 7o < Ts , a smg (Y,Ts), in the first term of eq. (4.14). 

The results for the rapidity spectra are shown in figures 7 and 9 for Drell-Yan and 
Higgs production, respectively. The two contributions from 7o < Ts and % > Ts are shown 
in red and green and the total result given by their sum in black. The error bars here show 
the scale variations up and down by a factor of two. Note that the relative size of the two 
contributions and the degree of cancellation between them can change significantly as the 
scale (or Ts value) is changed. To validate the results from 7o-slicing method, we compare 
to results from Vrap [106, 107] for Drell-Yan and from HNNLO [39, 110] for Higgs, which 
are shown by the blue line and band. For both processes we find excellent agreement. 

In figures 8 and 10 we show the fractional difference of the 7o-slicing results relative to 
Vrap and HNNLO, respectively. At NLO with Ts = 0.03 GeV, the agreement is excellent. 
For Drell-Yan at NNLO, there is a small offset between the two results visible in figure 8, 
representing 0.4% of the total cross section. A similar offset of —0.2% is also present 
in the Higgs case, but hardly visible because the scale variations are much larger. It is 
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Figure 7. The NNLO rapidity distribution in Drell-Yan production. We plot the various ingredients 
in the 7o-slicing method for Ts = 0.1 GeV, where in all cases the error bars correspond to the up 
and down scale variation. The blue histogram shows for comparison the NNLO result from Vrap. 
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Figure 8. The scale uncertainty band in the Drell-Yan rapidity distribution for both Vrap and 
7o-slicing, relative to the central scale from Vrap at NLO (right) and NNLO (left). 


due to the missing 0(<5ir) nonsingular terms for To < Ts = 0.1 GeV. A smaller value 
of Ts would be needed to reduce this effect. The size of the missing nonsingular terms 
we observe here is consistent with their expected size from our estimates in section 2.3.3. 
Nevertheless, it is actually encouraging to see that even with the simple 7o-slicing we are 
able to obtain this level of agreement. We would expect that an implementation of the 
differential 7o-subtractions will allow one to use <5 jr values well below 10~ 4 . 

We conclude this discussion by noting that it is important, particularly for more com¬ 
plex processes, to carefully quantify the size of the neglected 0(<5ir) nonsingular contribu¬ 
tions. In particular, as already seen in figure 1, one cannot draw any conclusions for their 
possible size at NNLO from knowing their size at NLO. Also, the difference in the result 
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Figure 9. The NNLO rapidity distribution in gg H production. We plot the various ingredients 
in the 7o-slicing method for Ts = 0.1 GeV, where in all cases the error bars correspond to the up and 
down scale variation. The blue histogram shows for comparison the NNLO result from HNNLO. 



Figure 10. The scale uncertainty band in the Higgs rapidity distribution for both HNNLO and 
7o-slicing, relative to the central scale from HNNLO at NLO (right) and NNLO (left). 


when varying the Ts value is not necessarily a good estimate of the absolute size of the 
missing nonsingular terms, because as discussed in section 2.3.3, their scaling with (5 ir for 
<5m —> 0 is much weaker than linear. A crucial check one should perform is to plot the 
nonsingular distribution as in figures 4 and 6 and check its convergence toward zero. 

5 Conclusions 

Higher-order computations in QCD require the use of some subtraction technique that 
allows one to extract the collinear and soft phase-space divergences from the real-emission 
diagrams, and cancel these against the explicit divergences from the virtual loop diagrams. 


- 35 - 










































We explained how a subtraction scheme can be constructed using an IR safe IV-jet resolution 
variable to control the approach to the IR-singular limit. The IV-jettiness observable Tn 
is ideally suited for this task due to its known and simple factorization properties. Our 
resulting IV-jettiness subtraction method is similar in spirit to the qx subtraction method 
introduced by Catani and Grazzini for color-singlet production, but may be applied to 
processes with arbitrarily many colored final-state partons (plus any color-singlet final 
state). 

In our method, the subtraction term corresponds to the appropriate fixed-order ex¬ 
pansion of the singular IV-jettiness cross section, which can be efficiently computed using 
SCET. In this context, SCET allows the subtraction term to be broken down into various 
pieces (beam, jet, and soft functions) that are easier to compute, with the beam and jet 
functions being reusable for processes with any number of jets. The extension to N 3 LO 
is possible and requires the calculation of the beam, jet, and soft functions at three-loop 
order. 

We discussed in depth the details of the subtraction procedure, giving explicitly the 
equations and ingredients needed to construct the 7jv-subtraction terms at NLO and 
NNLO. The only ingredient which is not explicitly known is the ^-independent constant 
term of the NNLO IV-jet soft function for three or more IV-jettiness axes. It can however 
be obtained relatively straightforward with existing technology. We also discussed how 
the IV-jettiness subtractions can be implemented in practice. To demonstrate the method 
and study some of its numerical aspects, we presented NNLO results for the Drell-Yan 
and Higgs rapidity spectra computed using O-jettiness subtractions in its simplest form 
as a slicing method. The slicing method has been previously shown to be successful for 
NNLO computations in ref. [69] and very recently in refs. [72, 73]. Given the viability of 
the 7}v-slicing, it will be very interesting to extend the implementations to the differential 
IV-jettiness subtractions. 

We have also suggested and discussed several different ways in which the numerical 
convergence of the IV-jettiness subtraction method can be systematically improved. One 
option would be to include the leading nonsingular terms in the subtraction. These correc¬ 
tions are described by subleading factorization theorems for IV-jettiness and SCET offers 
a systematic framework to compute them. Another way to improve the numerical con¬ 
vergence would be to make the subtraction more local, by performing the subtraction 
differentially in additional observables (such as px) and/or splitting the total IV-jettiness 
observable into its components in the jet and beam regions. Much of the recent work in 
SCET on deriving factorization formulae for multi-differential cross sections can be very 
useful in this direction. 

We only explicitly discussed the case of massless partons here. The construction of 
analogous 7Tv-subtractions for processes involving massive quarks is possible with the same 
techniques. For m q <C Q, one would consider a massive quark jet with its own IV-jettiness 
axis making use of the tools in SCET developed for the treatment of massive collinear 
quarks [114-119]. For Q ~ m q , e.g. tt pair production and similar processes, an analogous 
approach to refs. [48, 49] can be used. This amounts to treating the heavy quarks as part of 
the hard interaction (without its own IV-jettiness axis) together with a more complicated 
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soft function to account for soft gluon emissions from the heavy quarks. We leave further 
development in this direction to future work. 
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A Subtraction Ingredients 

We write the a s expansion of the QCD beta function and the cusp and noncusp anomalous 
dimensions as 

i 00 11 

(1 w—-v / r\ \ 77.+1 

/x—a s (/x) =/3[a s +)], /3(a s ) = -2a s (y J , (A.l) 

^ n =0 7F 

and 

r„, P K) = f;r„(g )" +1 , 7 M«,) = ET„(g)" +1 . ( a . 2 ) 

77—0 77—0 

The coefficients of the MS beta function and cusp anomalous dimensions we need are 

A) = y C A - g T f n f , 

Pi = y C\ - (y C A + 4 C F ) T f n f , (A.3) 

and 


r* = c F r n , 
r 0 = 4, 


= c A r n , 


r, =4 


Ca "'] = t [ CA(i - D + 5ft J 


For the quark jet and beam functions in MS we have [94, 96] 

7jo = 7bo = , 

7ji = 7b 1 = C F Ca(^~ - 8 OC 3 ) + C F { 3 — 4+ + 48£s) + /3o(y—I—y) 

For the gluon jet and beam functions in MS we have [95, 97, 120] 

7jo = 7bo = 2/3 0 5 

'94 2vr 2 — 




+ 2/3i. 


(A.4) 


(A.5) 


(A.6) 
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A.l Jet function 

We write the a s expansion of the quark (i = q) and gluon (i = g) jet functions as 


Ji(s,n) = Y 


n =0 


a s (v) 

47r 




s,g). 


The coefficients have the form 


2m—1 

4 m) (»,rt = .4”h(»)+ x 

^_n t 1 7 


n =0 


(A.7) 


(A.8) 


where the C n {x) are plus distributions as defined in eq. (2.8). The jet function is naturally 
a distribution in s/gA, and this is the only g dependence of the coefficients. Rescaling the 
arguments of the distributions using eq. (2.13), we have 


2 m —1 


r\Q,k,,) - t 41(f) *<«+f e 4:’(f) f-(f) ■ 


2m—1 

A m h\\ — A m ) \r( m ) 

~ 1 + J i,n n + i ’ 

n=0 

2m-l-n / 7 n. 

A m )f\\- r( m ) , ( n + r( m ) i,A 


Qi 
ln n+1 A 


c 0 (a) = 4:>+ e 


fc=i 


_ 

n j 7.1 j,n+fc 


In A. 


(A.9) 


where £ is an arbitrary dimension-one parameter, which exactly cancels between the differ¬ 
ent rescaled coefficients and that we can choose at our convenience. The J-'Y (X) are the 
coefficients appearing in the explicit expressions for the subtraction terms in sections 3.3.1 
and 3.3.2. 

The jet function coefficients in eq. (A.8) read up to two loops 


A 1 ) — p* 

J i, 1 — 1 0 > 
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(A.10) 
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The 5(s) pieces for the quark jet function are [121, 122] 


J (0) -1 

J q -1 — 1 > 

Jql-1 = C F (7 — IT 2 


Jq 2 ~l = C F 
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/205 
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and for the gluon jet function they are [95, 120, 123] 


j(°) — 1 

J a -1 “ 1 ’ 


J ^ i - CA (r ff2 ) + 3 A 


t(2) _ /~<2 

J g,~ 1 - L V1 


4 
3 

/4255 
V 108 
'25 


5 

26» 2 


+ 1^_ 72C3 ) +c , M _^_^ + ^) (a , 2) 


9 


180 

'55 


A.2 Beam function 

The beam function is given by [51, 96] 

Bi(t,x,n ) = ^ z, n, hf) fj (7, Mf) , 


(A.ll) 


115 65vr 2 56C 3 \ 


(A.13) 


where fj(x, /j,p) are the standard PDFs and Tjj(t, z, fi. /ip) are perturbative matching coef¬ 
ficients. Here, we have explicitly separated the fip dependence, which cancels between the 
matching coefficients and the PDFs, such that the beam function is [ip independent up to 
higher orders in a s (fi). (Usually, one takes fip = fi in the fixed-order beam function, since 
these are not really formally distinct scales.) For our purposes, the fip dependence in the 
beam function determines the complete fip factorization scale dependence in the singular 
fixed-order cross section, while the fi dependence contributes to the usual renormalization 
scale dependence. 

We expand the beam function matching coefficients as 

= ^2l^\t,z,n,n F ) (^~) • (A.14) 

77 .— 0 

The perturbative coefficients have the structure 

2m— 1 

m + E 

n=0 




z , —w I 
Hp J P 




(A.15) 


r (m) 




(' t,z,n,HP) =1^1^, 


4^ 

/4 
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where the C n {x) are the plus distributions defined in eq. (2.8). The beam function is nat¬ 
urally a distribution in f//r 2 . Rescaling the arguments of the distributions using eq. (2.13), 
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we have 
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4™>feA F ,A)=jN(2,A F )+ £ to*A, (A.16) 

fc=l 

where £ is an arbitrary dimension-one parameter, which exactly cancels between the differ¬ 
ent rescaled coefficients and that we can choose at our convenience. From these coefficients 
we also define the corresponding beam function coefficients as 

b'^W.w.a) = x/ *(? w ) ■ (AJ7) 

which are the coefficients appearing in the explicit expressions for the subtraction terms in 
sections 3.3.1 and 3.3.2. 

The results for the coefficients in eq. (A. 15) are as follows. At LO, we simply have 

4'-i (z,X F ) = S ij 6(l-z). 

The NLO coefficients have been computed in refs. [51, 96, 97], and are given by 


(A.18) 
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ij i 'vr 

The NNLO coefficients have been computed in refs. [59, 60], and read 
l%l(z,\F) = \(T i 0 ) 2 6 ij 5(l-z), 


(A.19) 
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(A.20) 


X^l l {z,\ F ) = u[f{z) + 
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Explicit results for the matching functions /^(z) and 1-^ (z) as well as the splitting func¬ 
tions P^\z) and P--p(z) and all required convolutions between them can be found in 
refs. [59, 60] in the same notation that we use here. 


A.3 Single-differential soft function 

The single-differential iV-jettiness soft function is related to the one of eq. (3.18), which is 
multi-differential in the soft contributions to the 7jy, by 


S K (k, {&}, n) = yj [PJ d h s(k - £ A;*) S K ({/n}, {<?*}, //). (A.21) 


Recall that the subscript k encodes the information on the Born partonic channel. For the 
soft function, it specifies the color space of the external partons in which it acts. 

We expand the soft function in a s (/i) as 

§.(*, {«},(*) = £(£)" {*}.«). (A. 22 ) 

n 

where the perturbative coefficients can be written as 
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(A.23) 


The soft function is naturally a distribution in k/fx and this is the only fx dependence of 
the coefficients. Rescaling the arguments of the plus distributions using eq. (2.13), we have 
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(A.24) 


The dimension-one parameter £ is again arbitrary and exactly cancels between the coeffi¬ 
cients. The coefficients ({%}, A) are those appearing in the explicit expressions for the 
subtraction terms in sections 3.3.1 and 3.3.2. In the rest of this subsection the dependence 
on the jet axes qi of the soft function and its anomalous dimension is always understood 
and we often suppress the explicit {<jj} argument. 
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The renormalization scale dependence of the soft function is subject to the renormal¬ 
ization group equation derived in ref. [52], 


M ) = \ [ d k' [%(fc - k') S K {k') + S K (k - k')i s {k') 
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with the soft anomalous dimension 
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+ 7s[u s (v)}d{k) 

= r cusp [a s (/i)]|2cE 0 (^) + S(k) [L({ljj}) + I] | + 7 sMm)] S(k ). (A.26) 

Here, A^ = 1 if the partons z and j are both incoming or both outgoing and A^ = 0 if 
one of them is incoming and the other one outgoing. The invariant 

Sij = = 2® • qj (A.27) 


is always positive with our conventions and corresponds to an angular measure between 
any two partons (depending on the precise choice of the Qi). Note that T cusp (a s ) here 
has the overall color factor removed, see eqs. (A.2) and (A.4). To write the last line in 
eq. (A.26), we defined the abbreviations 

C = E Tf = 1* E Ci ( with C <1 = C « = Cf ’ C 9 = Ca ) > 

i i 

L({%}) = E T ' : ' T i ln ®b' ’ 

I = ivE'I'. -T, A ij = ivr[2(T a + T b ) 2 - C] . (A.28) 

*7 


Note that for ee and ep collisions, I is always proportional to l re and can be ignored, as 
it drops out of eq. (A.25). Similarly, for pp collisions it can be ignored for 0-jet and 1-jet 
processes where the color space is still trivial. Up to two loops the noncusp soft anomalous 
dimension is given by 


7§{a«)=0 + C 7S i(g) 2 +O( Q »), 

7 S1 =a,(—f+ 28C,)-4). (A.29) 

The fixed-order coefficients in eq. (A.23) are as follows. At leading order, we have 

§i 0 h({«}) = l«. (A.30) 

The one-loop coefficients are given by [52] 


5$({«}) = - 2 r 0 c, 

^!o(te}) = -r 0 L({%}), 


§Ei(te}) = E T *- T ; 

i¥=j 



g +4 ^2 Iij,m({qi}) , 

m^i,j 


(A.31) 
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where 


hw.H*})-2e, f-2!-, «L,*4 . )ln: 

$ij &ij ^ Sjm Sim ' t z F' l i3i rn ' 


} jm 




?ij 


{: 


’j™ 


1 l^i,j,m) 


(A.32) 


The /o and I\ are finite phase-space integrals, which are required for three or more N- 
jettiness axes. They are not known fully analytically, but can be evaluated numerically 
for a given set {<)*}. Their explicit expressions and an algorithm to reduce them to simple 
one-dimensional numerical integrals for arbitrary N is provided in ref. [52]. With only 
three IV-jettiness axes, the integrals are still planar. Starting from four axes, the angles 
4>im also enter, which are the azimuthal angles between the q m and qi axes in the plane 
transverse to the (fi and qj axes. 

Iteratively solving the RGE in eq. (A.25), we obtain the two-loop coefficients 


££>({&}) = 2T 2 C 2 , 

S 2 2({®}) = roC[3r 0 L + 2/3 0 ], 

= rg(L 2 + x - [I, L] - ^C 2 ) + 2T 0 (A)L - C S™ !({&})) - 2^0 , 

2 

£$({%}) = rgc(4CCs - yL) - TiL - C 751 

- y ({L, £«,({%})} + [I, ^({ft})]) - 2ft ({ft}). (A.33) 


For two external partons, k = qq and n = gg , the result for the two-loop constant is 
known analytically [124-126] and does not depend on the whether the partons are incoming 
or outgoing, i.e., it is the same for 0 —> qq, q —> q, and qq —> 0, and similarly for two 
gluons [127], 


sgli = Cf 


s£!-i = c* 


Ct(- 


640 4tt 2 22vr 4 
’"27" + IT + “45~~ 


n ! 640 | 47t 2 177r 4 

+ “3” + “90" 


3vr 4 


/ 20 37vr 2 58Cs\ 

+/3°(^-y- — + —j 

„ , 20 37tt 2 

+ A)(_ 27~^ + 


58Ca\ 

3 ) 


(A.34) 
(A.35) 


The two-loop constants Sk (k = ggg, qqg) required e.g. for 1-jettiness in pp collisions, 
have recently been computed numerically in ref. [89]. The two-loop constant for arbitrary 
iV-jet processes can in principle be obtained numerically from known results for two-loop 
soft amplitudes as outlined in ref. [89]. 
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